import numpy as np
from scipy.integrate import quad
An integral that comes up when trying to find the chemical potential of a system of identical Bosons is:
$\int_0^\infty \dfrac{x^2}{e^{x^2-\eta}-1}\, dx$
where $\eta=\mu/k_\mathrm{B}T$ and $\mu$ is the chemical potential.
This short tutorial deminstrate how to use quad() to numberically evaluate this integral for the $\eta = 0$ case.
def integrand(x):
return x**2/(np.exp(x**2) - 1)
y, err = quad(integrand, 0, np.inf)
<ipython-input-27-472f14518a59>:2: RuntimeWarning: overflow encountered in exp return x**2/(np.exp(x**2) - 1)
y, err
(1.1575786866969464, 8.140518531748718e-09)
y, err = quad(integrand, 0, 10)
y, err
(1.1575786866970472, 7.989498536842354e-12)
For this choise of the upper integration limit, there was no warning and our result for y changed by only a completely negligible amount.